
\chapter{数值求导、多项式拟合}

\section{数值求导}
我们常常需要计算一个函数的导数，
数值的做法（永远）是寻求近似。
\subsection{导数的数值近似}
于是我们写出函数的泰勒展开
\begin{eqnarray}
f(a+h) &=& f(a) + hf'(a) + \frac{f''(a)}{2}h^2 + \frac{f'''(a)}{6} h^3 +\cdots \\
f(a-h)  &=& f(a) - hf'(a) + \frac{f''(a)}{2}h^2 - \frac{f'''(a)}{6} h^3 +\cdots
\end{eqnarray}
其中，$h$ 是一个很小的数。
我们可以用如下公式计算 $f'(a)$，
\begin{eqnarray}
	f'(a) = \frac{ f(a+h) - f(a) }{h} + O(h), \\
	f'(a) = \frac{ f(a) - f(a-h) }{h} + O(h).
\end{eqnarray}
所以，我们计算
\begin{eqnarray}
	f'(a) \approx \frac{ f(a+h) - f(a) }{h} , \\
	f'(a) \approx \frac{ f(a) - f(a-h) }{h} .
\end{eqnarray}
作为 $f'(a)$ 的近似，误差为 $O(h)$。

我们也可以设计
\begin{eqnarray}
 f'(a) = \frac{ f(a+h) - f(a-h) }{2h} + O(h^2),
\end{eqnarray}
并用相应的近似公式来计算导数
\begin{eqnarray}
f'(a) \approx \frac{ f(a+h) - f(a-h) }{2h} , 
\label{formula:derivative}
\end{eqnarray}
这个公式的误差就是 $O(h^2)$，所以这个公式更精确。

如果我们需要两阶导数，我们也可以类似地构造
\begin{eqnarray}
	f''(a) = \frac{ f(a+h) + f(a-h) - 2 f(a) }{h^2} + \frac{ f^{(4)} (a) }{12} h^2 + O(h^4),
	\label{double-derivative}
\end{eqnarray}
用近似公式来计算 $f''(a)$，
\begin{eqnarray}
f''(a) \approx \frac{ f(a+h) + f(a-h) - 2 f(a) }{h^2},
\end{eqnarray}
其误差约为 $ \frac{ f^{(4)} (a) }{12} h^2 $。

一个自然的想法是：我让 $h$ 特别小，$h^2$ 也会特别小，那么我的导数岂不是特别精确了吗？
然而并不是这样的。
当 $h$ 特别小的时候，$f(a+h)$ 与 $f(a-h)$ 会特别接近，就出现了前面提到的“两个非常相近的数相减”的问题。

下面我们写一个程序，计算 $e^x$ 在 $x=10$ 处的二阶导数，来观察舍入误差的影响。
我们定义相对误差为
\begin{equation}
	\epsilon_r (x)= \frac{ f''_{cal} (x) - f''_{real} (x) } { f''_{real} (x)},
\end{equation}
其中 $f''_{cal} (x), f''_{real} (x)$ 分别为二阶导数的计算值、真实值。
\begin{lstlisting}
#include<iostream>
using namespace std;

#include<cmath>//包括cmath库，后文中就可以调用库中的exp()函数

int main(){

	double h[10];//步长
	for(int i=0;i<10;i++){
		h[i] = pow(0.1,i+1);
        //　10个步长分别为0.1,0.01,...,1E-10
	}

	double x = 10;
	double double_derivative;//两阶导数 f''(x) 的计算值
	double relative_error;//计算值的相对误差

	FILE *fp=fopen("double_derivative_ex.txt","w");//定义一个文件流，打开一个 txt 文件，后面利用这个 txt 文件作图
	if(fp==NULL){//如果打开文件失败
		cout<<" failed to open double_derivative_ex.txt"<<endl;
		exit(1);//结束程序，返回异常信号，exit(1)表示异常退出，exit(0)表示正常退出
	}

	for(int i=0;i<10;i++){
		double_derivative = ( exp(x+h[i]) + exp(x-h[i]) - 2*exp(x) )/h[i]/h[i]; //计算两阶导数
		relative_error = double_derivative/exp(x) -1 ; //计算相对误差
		cout<<" h= "<<h[i]<<"  relative error = "<< relative_error <<endl; //输出到屏幕上
		fprintf(fp, "%lf \t %lf \n", log10( h[i] ), log10( fabs(relative_error) ) ); //输出到文件里,用于画图　\t表示一个tab的空白，\n表示换行
	}
	fclose(fp);//关闭文件

	return 0;
}
\end{lstlisting}
在上面的程序中，我们用了c++中的FILE类型，即文件流，fp是一个文件流指针。
fopen　函数打开一个文件，"w"表是写入内容，"r"表示读取内容。
如果 fopen　打开文件失败，会返回NULL值，否则返还文件流指针，通过这个指针，就可以用fscanf或fprintf函数读取或写入内容。

我们将 $ \log_{10} | \epsilon_r | ~ \log_{10}(h)$ 做成图片，可以直观地看到相对误差与步长的关系。

\begin{figure}
	\centering
    \input{double_derivative_ex.tex}
    %\includegraphics[width=0.8\textwidth]{double_derivative_ex.eps}
	\caption{$ \log_{10} | \epsilon_r | \sim \log_{10}(h)$： $e^x$ 在 $x=10$ 处的两阶导数的相对误差与步长的关系}
	\label{double-derivative-relative-error}
\end{figure}
我们发现，相对误差在 $h \approx 10^{-4}$ 时，取得最小值。
怎样理解这个现象？

\subsection{误差分析}
误差分为两部分，一部分是截断误差，一部分是舍入误差。
\begin{equation}
	\epsilon = \epsilon_{truncation} + \epsilon_{round-off},
\end{equation}
截断误差即公式\ref{double-derivative}中被舍去的部分，大约为
\begin{equation}
	\epsilon_{truncation} \approx \frac{ f^{(4)} (a) }{12} h^2,
\end{equation}
而舍入误差，按我们前面“数值误差”一章的分析，
\begin{equation}
	\epsilon_{round-off} \approx \frac{ 2f(a) }{ h^2} 10^{-16},
\end{equation}
所以总误差为
\begin{equation}
	\epsilon = \frac{ f^{(4)} (a) }{12} h^2 + \frac{ 2 f(a) }{ h^2} 10^{-16},
\end{equation}
将它对 $h$ 求导数，导数为零时误差取得最小值，得到 $ d \epsilon / dh = 0$ 时，
\begin{equation}
	h \approx 2.2 \times 10^{-4},
\end{equation}
所以 $h$ 取 $2.2 \times 10^{-4}$ 时，计算误差最小。
这与我们在图 \ref{double-derivative-relative-error} 中看到的现象是一致的！

所以，我们想通过不断减小$h$的值，来得到尽量精确的导数，这个策略是行不通的。

于是人们开发了其他的策略，比如Richardson外推法。

\section{Richardson外推法求导}
前面我们追求数值求导的精度时，遇到了困难。
我们发现，步长太小时，舍入误差太大了，计算可能失去意义。
这一节，我们介绍Richardson的外推法，这是一个数值求导算法。

我们注意到，用(\ref{formula:derivative})求得的数值导数有如下形式，
\begin{equation}
D(h) = D(0) + a_2 h^2 + a_4 h^4 + \cdots,
\end{equation}
其中$D(h)$表示步长选为h时求得的数值导数，而$D(0)$表示导数的理论值。
常数$a_2,a_4,\cdots$与$h$并没有关系。

那么，我们可以定义
\begin{eqnarray}
D^0_0 &=& D(h) = D(0) + a_2 h^2 + a_4 h^4 + \cdots,\\
D^1_0 &=& D(h/2) = D(0) + a_2 (h/2)^2 + a_4 (h/2)^4 + \cdots,\\
&\cdots& \nonumber\\
D^k_0 &=& D(h/2^k) = D(0) + a_2 (h/2^k)^2 + a_4 (h/2^k)^4 + \cdots.
\end{eqnarray}
所以你看得出来，利用上面的各行，消掉$a_2, a_4, \cdots$就可以得到越来越精确的导数。
我们可以定义
\begin{eqnarray}
D^k_m = D^{k+1}_{m-1} + \frac{ D^{k+1}_{m-1} - D^k_{m-1} }{ 4^m -1 }, 
\end{eqnarray}
可以利用这个公式，递推出下表中的D值\\
\begin{center}
\begin{tabular}{cccc}
 $D^0_0$ 	&  &  &  \\  
 $D^1_0$	&	$D^0_1$ &  &  \\  
 $D^2_0$	&	$D^1_1$ & $D^0_2$  	&  \\  
 $D^3_0$	&  	$D^2_1$	& $D^1_2$ 	& $D^0_3$  \\  
 $\cdots$	&	$\cdots$  & $\cdots$ 	& $\cdots$ \\ 
\end{tabular} 
\end{center}
上图中，第一列具有$O(h^2)$精度，第二列具有$O(h^4)$精度，第三列具有$O(h^6)$精度。
这可以观察出来，我们仅展示$D^0_1$
\begin{eqnarray}
D^0_1 &=& D^1_0 + \frac{ D^1_0 - D^0_0 }{ 4^1 -1 } \nonumber \\
      &=& D(0) - \frac{1}{4} a_4 h^4 - \frac{5}{16} a_6 h^6 + \cdots.
\end{eqnarray}

写成代码如下
\begin{lstlisting}
/*
 * Richardon's algorithm to extrapolate an accurate derivative of func at x
 * It calculates naive derivatives using formula
 *      D^k_0(h/2^k) = (f(x+h/2^k) - f(x-h/2^k))/ ( h/2^{k-1} ), k=0,1,...,n
 * and then extrapolate an accurate derivative and returns it.
 * 
 * n mustn't be less than 0, otherwise it refuses the task and exit
 * 
 * Error:        O(h^{2(n+1)}
 */
double Richardson_deriv(double x, double h, int n, double (*func)(double) ){
        
        double * deriv = new double [ n+1 ];
        double coef,y,a,b,c;
        for(int i=0;i<=n;i++){
                coef = pow(0.5,i);
                deriv[i] = ( func(x + h*coef ) - func(x - h*coef) ) / (2* h* coef);
        }       
        for(int i=1;i<=n;i++){
                a = deriv[i-1];
                b = deriv[i];
                for(int j=i;j<=n;j++){
                        c = b + (b-a)/( pow(4.0,i) -1 );
                        a=b;
                        if(j<n)b=deriv[j+1];
                        deriv[j] = c;
                }
        }
        //cout<<" error is O(h^"<<2*n+2<<")\n";
        y = deriv[n];
        delete [] deriv;
        return y;
}
\end{lstlisting}
用这个办法，求取$exp(x)$在$x=0,1,2,3,4,5$处的一阶导数，得到的结果如表\ref{table:Richardson}，可以看到，仅仅推到$D^0_2$，误差已经小到$10^{-10}$量级。
\begin{table}[h]
\centering
\begin{tabular}{cccccc}
\hline 
$x$ & $D^0_0$ & $D^0_1$ & $D^0_2$ & error of $D^0_2$  \\ 
\hline 
0 & 1.00166750 & 0.99999979 & 1.00000000 & 0.00000000000 \\ 
1 & 2.72281456 & 2.71828126 & 2.71828183 & 0.00000000001   \\ 
2 & 7.40137735 & 7.38905456 & 7.38905610 & 0.00000000002   \\ 
3 & 20.11902956 & 20.08553274 & 20.08553692 & 0.00000000006   \\ 
4 & 54.68919246 & 54.59813866 & 54.59815003 & 0.00000000017   \\ 
5 & 148.66063807  & 148.41312817 & 148.41315910 &  0.00000000046  \\ 
\hline 
\end{tabular} 
\label{table:Richardson}
\caption{数值求导，Richardson算法：exp(x)的一阶导数}
\end{table}

\section{多项式插值}
如果我们已经知道了$N+1$个数据点$y_0 = f(x_0), y_1 = f(x_1), \cdots, y_N = f(x_N)$，其中$x_i$都不相同。
但我们不知道函数$f(x)$的形式，却需要推测$f(x_{N+1})$的值，$x_{N+1}$有可能在两个$x_i$之间（即插值），也有可能在$x_i$覆盖的范围之外（即外推）。

我们可以用一个多项式，去拟合这$N+1$个已知的数据点，然后根据这个多项式，去做插值和外推。
多项式的形式为
\begin{equation}
P(x) = a_0 + a_1 x + \cdots + a_N x^N,
\end{equation}
所以我们需要用已知的数据点，来求参数$a_0, a_1, \cdots, a_N$。

你当然可以将$x_i$代进去，命令
\begin{equation}
P(x_i) = y_i,
\end{equation}
得到$N+1$个线性方程，求解这个线性方程组，来完成任务。

\paragraph{拉格朗日插值公式}
好在我们不需要费那个劲，拉格朗日已经给出了结果。
拉格朗日公式为
\begin{equation}
P(x) = \sum^N_{i=0} y_i \prod_{k \neq i} \frac{ x-x_k }{ x_i - x_k },
\end{equation}
很容易看出，$P(x)$经过$N+1$个已知数据点，而且是$N$阶多项式，所以正好满足要求！
因为我们知道，经过$N+1$个数据点的$N$阶多项式是唯一的。

误差可以计算（unchecked）
\begin{equation}
f(x) = P(x) + \frac{1}{(N+1)!} (x-x_0)\cdots(x-x_N)f^{(N+1)}(\epsilon),
\label{polynomial-error}
\end{equation}
其中$\epsilon$为包括$x_0,\cdots,x_N,x$的最小区间内的一个数。

\paragraph{牛顿插值公式}
牛顿也有一个递推公式，用来计算$P(x)$。

如果
\begin{equation}
P_{N-1}(x) = a_0 + a_1 x + \cdots + a_{N-1} x^{N-1}
\end{equation}
通过
\begin{equation}
(x_0, y_0), (x_1, y_1), \cdots, (x_{N-1}, y_{N-1}),
\end{equation}
N个点，则有
\begin{eqnarray}
P_N(x) &=& P_{N-1}(x) + a_N x^{N}, \nonumber\\
a_N &=& \left\{
\begin{aligned}
&0 &if ~~x_N=0,\\
&\frac{ y_N - P_{N-1}(x_N) }{a_N} &else
\end{aligned}
\right.
\end{eqnarray}
$P_N(x)$通过
\begin{equation}
(x_0, y_0), (x_1, y_1), \cdots, (x_{N-1}, y_{N-1}),　(x_N, y_N)
\end{equation}
N+1个点。
这是一个递推公式，逐渐增加数据点的个数，多项式逐渐增加高阶项，原来的低阶项并不变。

\section{Nevill外推法}

如果给定$x$的值，用拉格朗日公式，需要$N^3$次操作（加减乘除），不是一个很高效的算法。
所以Nevill给了一个更高效的算法，只需要$N^2$次操作。
我们定义
\begin{equation}
P_0, P_1, P_2, \cdots, P_N,
\end{equation}
为$N+1$个０阶多项式，分别通过第$0,1,\cdots,N$个点。
很显然，０阶多项式就是常数，
\begin{equation}
P_0 = y_0, P_1 = y_1, \cdots, P_N = y_N,
\end{equation}
而
\begin{equation}
P_{01}, P_{12}, P_{23}, \cdots, P_{N-1,N},
\end{equation}
为$N+1$个1阶多项式，分别通过第${0,1},{1,2},{2,3}\cdots,{N-1,N}$个点。
相应地
\begin{equation}
P_{i(i+1)\cdots(i+m)}
\end{equation}
为m阶多项式，分别通过第${i,i+1,i+m}$个点。

显然，我们真正需要的是$P=P_{01\cdots N}$。
然而，这些多项式具有递推关系
\begin{equation}
P_{i(i+1)\cdots(i+m)} = \frac{ (x-x_{i+m}) P_{i(i+1)\cdots(i+m-1)} - (x_i -x) P_{(i+1)\cdots(i+m) } }{ x_i - x_{i+m} },
\end{equation}
很容易证明等式右边确实通过第${i,i+1,i+m}$个点，我省去这个证明，请读者你自己练习。

举$N=3$为例，我们按下图推得$P_{0123}$，

\begin{tabular}{cccc}
$P_0$	&  &  &  \\   
$P_1$	&$P_{01}$&  &  \\   
$P_2$	&$P_{12}$&$P_{012}$ &  \\  
$P_3$	&$P_{23}$&$P_{123}$ & $P_{0123}$ \\ 
\end{tabular} 

给定$N+1$个数据点，给定$x$值，按如下步骤计算$P(x)$，
\begin{itemize}
    \item	计算$P_0(x),\cdots,P_3(x)$的值。
    \item	用$P_0,P_1$计算$P_{01}(x)$的值，类似地，计算第二列所有多项式的值。
    \item	计算第三列。
    \item 	计算第四列，即得到$P(x)$的值。
\end{itemize}

于是我们可以写一个代码，实现Nevill算法
\begin{lstlisting}
/*
 * Nevill algorithm
 *
 * Given n data pointes, in (x_data, y_data), and an arbitrary x value, a polynomial P(x) passes these n points, the function Nevill returns y=P(x)
 *
 * If the datapoints' x values are too close, round off error can be significant, therefore the function refuses the task
 */
double Nevill(int n, double *x_data, double *y_data, double x){
        double y;
        //check if data points are too close
        for(int i=0;i<n;i++){
                for(int j=i+1;j<n;j++){
                        if( fabs(x_data[i]-x_data[j]) < 1E-6 ){
                                cout<<" datapoints' x values are too close.\n";
                                exit(1);
                        }
                }
        }

        double * tempP = new double [n];
        for(int i=0;i<n;i++) tempP[i] = y_data[i];// initially P_0, P_1, \cdots, P_{n-1}
        //use the recurrence relations to derive the final polynomial
        for(int i=1;i<n;i++){
                for(int j=i;j<n;j++){
                        tempP[j] = (x-x_data[j])*tempP[j-1] + (x_data[j-i] - x)*tempP[j];
                        tempP[j] /= x_data[j-i] - x_data[j];
                }
        }

        y = tempP[n-1];
        delete [] tempP;//删除动态内存

        return y;
}
\end{lstlisting}
我们可以写一个主函数，调用上面写的Nevill函数，实现如下功能：
\begin{itemize}
\item 程序自动从一个数据文件中读取数据
\item 程序自动进行多项式拟合
\item 程序自动给出$x_0,\cdots,x_n$区域内的函数值，并输出到文件中，供画图用。
\end{itemize}
代码如下
\begin{lstlisting}
int main(){

        int n;//数据个数
        double x,y,dy;

        double *x_data;
        double *y_data;

        FILE *fp=fopen("Nevill_input.txt","r");
        if(fp==NULL){
                cout<<"error: failed to open Nevill_input.txt.\n";
                exit(1);
        }

        //读一遍输入文件，并将数据行数输出到屏幕上
        n=0;
        while( fscanf(fp, "%*c%*[^\n]") != EOF ){
                n++;
        }
        cout<<"n="<<n<<endl;

        //开辟动态内存，用于存放数据
        x_data = new double [n];
        y_data = new double [n];
        //读取数据
        fp=fopen("Nevill_input.txt","r");
        for(int i=0;i<n;i++){
                fscanf(fp, "%lf %lf", &x_data[i], &y_data[i]);
        }
        fclose(fp);//结束读取

        //为了画图，选取x轴范围
        double xmin = x_data[0], xmax = x_data[0];
        for(int i=0;i<n;i++){
                if( xmin > x_data[i] ) xmin = x_data[i];
                if( xmax < x_data[i] ) xmax = x_data[i];
        }
        xmin -= 1;
        xmax += 1;

        int num_grid = 100;//曲线上的点个数
        double step = (xmax - xmin)/num_grid;//曲线上的点的间距

        //计算每个x点的函数值，并输出到文件中
        fp=fopen("Nevill_curve.txt","w");
        if(fp==NULL){
                cout<<"error: failed to open Nevill_curve.txt.\n";
                exit(1);
        }
        for(int i=0;i<num_grid;i++){
                x = xmin + step*i;
                Nevill(n,x_data,y_data,x,&y);
                fprintf(fp,"%lf   %lf\n", x, y);
        }
        fclose(fp);//结束输出

        return 0;
}
\end{lstlisting}
当然，在main函数前面需要包括输入输出库，声明标准库命名空间，以及声明上面定义的Nevill函数，我们省去了这些内容。

如果我们把数据文件"Nevill\_input.txt"中的内容设为
\begin{lstlisting}
1	2
2	3
3	5
4	2
\end{lstlisting}
然后运行上面的程序，得到数据文件以后，可以用gnuplot作图，得到图\ref{figure:Nevill}，可以看到，拟合出来的曲线通过每一个点。
利用这条曲线，就可以进行插值和外推。
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{Nevill}
\label{figure:Nevill}
\caption{多项式拟合与插值：内维尔算法}
\end{figure}

也可以稍简化算法，定义
\begin{eqnarray}
C_{m,i} \equiv P_{i\cdots(i+m)} - P_{i\cdots(i+m-1)}, \\
D_{m,i} \equiv P_{i\cdots(i+m)} - P_{(i+1)\cdots(i+m)},
\end{eqnarray}
容易推得
\begin{eqnarray}
D_{m+1,i} = \frac{ (x_{i+m+1} - x) ( C_{m,i+1} - D_{m,i}) }{ x_i - x_{i+m+1}}, \\
C_{m+1,i} = \frac{ (x_i -x) ( C_{m,i+1} - D_{m,i} ) }{ x_i - x_{i+m+1} }
\end{eqnarray}
例如四个点的多项式插值，可以由下图逐步推得

\begin{table}[h]
\centering
\begin{tabular}{cccc}
$P_0$	&  &  &  \\  
        &$P_{01}$&  &  \\    
$P_1$	&        &$P_{012}$  &  \\   
    	&$P_{12}$&          & $P_{0123}$   \\  
$P_2$	&         &$P_{123}$ &  \\ 
        &$P_{23}$&          &  \\
$P_3$   &       &       &        
\end{tabular} 
\end{table}

代码如下
\begin{lstlisting}
/*
 * Nevill 算法
 *
 * 给定 n 个数据点　(x_data[], y_data[])，指定内插/外推点位置　x，Nevill()函数自动进行多项式拟合，得到　x　处插值结果 y，并计算误差 dy
 *
 * 如果数据点的　x　坐标太靠近，误差可能会比较大，所以程序中断。
 *
 */
double Nevill(int n, double * x_data, double * y_data, double x, double &y, double & dy){

        double *c = new double [n];
        double *d = new double [n];

        int i, m, ns=0;
        double dis=fabs(x-x_data[0]), dif, den;
        for(int i=0;i<n;i++){
                if( dif=fabs(x-x_data[i]) < dis ){// 找到离x最近的数据点x[ns]，从它的函数值开始逐步推测y
                        dis = dif;
                        ns = i;
                }
                c[i] = y_data[i];//赋初值　c,d,　可以证明，这个初值带入公式，确然得到第二列　P_{01}, P_{12}, ...
                d[i] = y_data[i];
        }
        y = y_data[ns];//y的初始值

        for(m=1;m<n;m++){
                for(i=0;i<n-m;i++){
                        den = x_data[i] - x_data[i+m];
                        if( fabs(den) < 1E-9 ){
                                cout<<" two data points' x coordinates are too close.\n";
                                exit(1);
                        }
                        den = ( c[i+1] - d[i] )/ den;

                        d[i] = (x_data[i+m] - x) * den;//迭代
                        c[i] = (x_data[i] - x) * den;
                }
                y += ( dy = (2*ns < n-m ? c[ns]: d[--ns]) );// 这一行保证推测路径最接近直线。 ns < (n-m)/2，说明 ns 在上半区，由　c路径 （斜向下），否则由　d路径（斜向上）
        }
        delete [] c;
        delete [] d;
}
\end{lstlisting}
